_ UNCLASSIFIED 

SECuniTv  ci.  ASS'^'C  »TiON  Or  T*"S  °»C 

reporFdocumI 

'  REPORT  NUMBER  ”” 

AI  Memo  1114 

4  TITLE  find  Submit) 


?EAD  INSTRUCTIONS 

AD-A212  457 


Parallel  and  Deterministic  Algorithms  from  MRFs: 
Surface  Reconstruction  and  Integration 


7.  AUTHOR!.; 

Davi  Geiger  and  Federico  Girosi 


-•  REPORT  »  PERIOD  COVERED 

memorandum 

«.  "performing  org.  report  numier 


I.  CONTRACT  OR  GRANT  NUUIERflJ 

DACA76-85-C-0010 

N00014-85-K-0124 


tO.  PROGRAM  ELEMENT,  PROJECT.  TASK 
AREA  •  WORK  UNIT  NUMBERS 


.  PERFORMING  ORGANIZATION  NAME  AND  AOORESS 

Artificial  Intelligence  Laboratory 
545  Technology  Square 
Cambridge,  MA  02139 


II.  CONTROLLING  OPPICE  NAME  ANO  AOORESS 

Advanced  Research  Projects  Agency 
1400  Wilson  Blvd. 

Arlington,  VA  22209  _ 


14  MONITORING  AGENCY  NAME  A  AOORESSflf  dlllmrmtl  Iran  Canfrallln*  OHIca;  II.  SECURITY  CLASS,  fal  ihla  rapart; 

Office  of  Naval  Research  UNCLASSIFIED 

Information  Systems 
Arlington,  VA  22217 


IZ.  REPORT  DATE 

May  1989 


IS.  NUMBER  OP  PAGES 

37 


II*.  OECL  A  SSI  PIC  ATI  ON  /  DOWN  G  R  AOING 
SCHEDULE 


II.  DISTRIBUTION  STATEMENT  fal  (Ala  ftapari; 

Distribution  is  unlimited 


DTiC 

ELECTE 
AUG  1 5  1989 


'Tfcg, 


is.  supplementary  notes 

None 

IS.  KEY  WORDS  fCanllmia  an  raaaraa  alaa  II  na< 

eaaaarr  aN  IWanill?  A,  SJacA  maaSarJ 

surface  reconstruction 

parameter  estimation 

Markov  random  fields 
mean  field 

'deterministic  algorithms 

integration 

10.  ABSTRACT  fCMUfmM  «i  If  m4  ItfOTtlff  Of  M«c* 

In  recent  years  many  researchers  have  investigated  the  use  of  Markov  random  fields 
(MRFs)  for  computer  vision.  They  can  be  applied  for  example  in  the  output  of  the  vi¬ 
sual  processes  to  reconstruct  surfaces  from  sparse  and  noisy  depth  data,  or  to  integrate 
early  vision  processes  to  label  physical  discontinuities.  Drawbacks  of  MRFs  models 
have  been  the  computational  complexity  of_the  implementation  and  the  difficulty  in 
estimating  the  parameters  of  the  model.  ...  / 

I  In  this  paper  we  derivd*deterministic  approximations  to  MRFs  models.  One  of  the 


DD  ,jan“z  1473  EDITION  OP  I  NOV  III!  OBSOLETE  UNCLASSIFIED 

S/N  o:oi-om*  1101  I  -  - - 

SECURITY  CL  ASSIPICATION  OP  this  PAGE  f*»»an  Data  Brim 


'^o.  «'vuuX 

considered  models  is  shown  to  give  in  a  natural  way  the  graduate  non  convexity  (GNC) 
algorithm, proposed  by  Blake  and  Zisserman.  ^This  model  can  be  applied  to  smooth 
a  field  preserving  its  discontinuities.  A  new  model  is  then  proposed:  it  allows  the 
gradient  of  the  field  to  be  enhanced  at  the  discontinuities  and  smoothed  elsewhere.  All 
the  theoretical  results  are  obtained  in  the  framework  of  the  mean  field  theory,  that  is 
a  well  known  statistical  mechanics  technique.  A  fast,  parallel  and  iterative  algorithm 
to  solve  the  deterministic  equations  of  the  two  models  is  presented,  together  with 
experiments  on  synthetic  and  real  images.  The  algorithm  is  applied  to  the  problem  of 
surface  reconstruction  is  in  the  case  of  sparse  data.  We~ai50"~descri£e  a  fast  algorithm 
that  solves  the  problem  of  aligning  the  discontinuities  of  different  visual  models  with 
intensity  edges  via  integration.  /  ^ — 


j~Acce>io;  For 

NTlS  CRA&I 
DTlC  TAB 

Unannouuced 
JjStifiCdtton  . 

By _ _ _ 

Distribution! 


□ 

□ 


Availability  Codes 


Dist 


Avail  o  -d  |  of 
Special 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
ARTIFICIAL  INTELLIGENCE  LABORATORY 

and 

CENTER  FOR  BIOLOGICAL  INFORMATION  PROCESSING 
WHITAKER  COLLEGE 

A. I.  Memo  No. 1114  June  1989 

C.B.I.P.  Paper  No.  35 

Parallel  and  deterministic  algorithms  for  MRFs: 
surface  reconstruction  and  integration 

Davi  Geiger  and  Federico  Girosi 
Abstract 

In  recent  years  many  researchers  have  investigated  the  use  of  Markov  random  fields 
(MRFs)  for  computer  vision.  They  can  b-  applied  for  example  in  the  output  of  the  vi¬ 
sual  processes  to  reconstruct  surfaces  from  sparse  and  noisy  depth  data,  or  to  integrate 
early  vision  processes  to  label  physical  discontinuities.  Drawbacks  of  MRFs  models 
have  been  the  computational  complexity  of  the  implementation  and  the  difficulty  in 
estimating  the  parameters  of  the  model. 

In  this  paper  we  derive  deterministic  approximations  to  MRFs  models.  One  of  the 
considered  models  is  shown  to  give  in  a  natural  way  the  graduate  non  convexity  (GNC) 
algorithm  proposed  by  Blake  and  Zisserman.  This  model  can  be  applied  to  smooth 
a  field  preserving  its  discontinuities.  A  new  model  is  then  proposed:  it  allows  the 
gradient  of  the  field  to  be  enhanced  at  the  discontinuities  and  smoothed  elsewhere.  All 
the  theoretical  results  are  obtained  in  the  framework  of  the  mean  field  theory,  that  is 
a  well  known  statistical  mechanics  technique.  A  fast,  parallel  and  iterative  algorithm 
to  solve  the  deterministic  equations  of  the  two  models  is  presented,  together  with 
experiments  on  synthetic  and  real  images.  The  algorithm  is  applied  to  the  problem  of 
surface  reconstruction  is  in  the  case  of  sparse  data.  We  also  describe  a  fast  algorithm 
that  solves  the  problem  of  aligning  the  discontinuities  of  different  visual  models  with 
intensity  edges  via  integration. 

©  Massachusetts  Institute  of  Technology,  1989 
This  paper  describes  research  done  within  the  Center  for  Biological  Information  Processing, 
in  the  Department  of  Brain  and  Cognitive  Sciences,  and  at  the  Artificial  Intelligence  Labo¬ 
ratory.  This  research  is  sponsored  by  a  grant  from  the  Office  of  Naval  Research  (ONR),  Cog¬ 
nitive  and  Neural  Sciences  Division;  by  the  Alfred  P.  Sloan  Foundation;  by  the  Artificial  In¬ 
telligence  Center  of  Hughes  Aircraft  Corporation  (S 1-80 1 534-2 ) ;  and  by  the  NATO  Scientific 
Affairs  Division  (0403/87).  Support  for  the  A.  I.  Laboiatory’s  artificial  intelligence  research 
is  provided  by  the  Advanced  Research  Projects  Agency  of  the  Department  of  Defense  un¬ 
der  Army  contract  DACA76-85-C-0010,  and  in  part  by  ONR  contract  N00014  85-K  0124. 


O  ^  c.  Q  £ 


1  Introduction 


In  order  to  give  a  viewer  information  about  a  three  dimensional  scene  many 
algorithms  have  been  developed  on  several  early  vision  processes,  such  as  edge 
detection,  stereopsis,  motion,  texture,  and  color  .  This  information  refers  to 
properties  of  the  scene  as  shape,  distance,  color,  shade  or  motion,  and  it  is 
usually  noisy  and  sparse:  more  processing  is  then  necessary  to  extract  the 
relevant  information  and  fill  in  sparse  data.  In  recent  years  many  researchers 
[9][12][6]  [4] [5]  have  investigated  the  use  of  Markov  Random  Fields  (MRFs) 
for  early  vision.  MRFs  models  can  generally  be  used  for  the  recon'^rv.ction 
of  a  function  starting  from  a  set  of  noisy  sparse  data,  such  as  intensity, 
stereo,  or  motion  data.  They  have  also  been  used  to  integrate  early  vision 
processes  to  label  physical  discontinuities.  Two  fields  are  usually  required 
in  the  MRFs  formulation  of  a  problem:  one  represents  the  function  that 
has  to  be  reconstructed,  and  the  other  is  associated  to  its  discontinuities. 
The  essence  of  the  MRFs  model  is  that  the  probability  distribution  of  the 
configuration  of  the  fields,  given  a  set  of  data,  is  given  as  a  Gibbs  distribution. 
The  model  is  then  specified  by  an  “energy  function”,  that  can  be  modeled  to 
embody  the  a  priori  information  about  the  system.  In  the  standard  approach 
an  estimate  of  the  field  and  its  discontinuities  is  given  by  the  configuration 
that  maximizes  the  probability  distribution,  or  equivalently  that  minimizes 
the  energy  function.  Since  the  discontinuity  field  is  a  discrete  valued  field  (it 
assumes  only  the  values  0  or  1)  this  becomes  a  combinatorial  optimization 
problem,  that  can  be  solved  by  methods  of  the  Monte  Carlo  type  (simulated 
annealing[10],  for  example). 

The  MRFs  formulation  is  appealing  bee  use  describes  a  system  by  local  in¬ 
teractions  and  allows  to  capture  many  features  of  the  system  of  interest  by 
simply  adding  appropriate  te  ,  ,s  m  the  energy  function.  However,  it  has  two 
main  drawbacks:  the  amount  nputer  time  needed  for  the  implementa¬ 

tion  and  the  difficulty  in  estimating  the  parameters  that  control  the  relative 
weight  of  the  various  terms  of  the  energy  function. 

In  this  paper  we  propose  a  deterministic  approach  to  MRFs  models.  It 
consists  in  explicitly  writing  down  a  set  of  equations  from  which  we  can 
compute  estimates  of  the  mean  values  of  the  field  /  and  the  line  process. 
We  study  two  MRFs  models,  the  second  being  an  extension  of  the  first,  and 
we  show  that  in  each  case  the  deterministic  equations  lead  to  a  fast,  parallel 
and  iterative  algorithm  that  can  be  used  to  obtain  estimates  of  the  fields 
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everywhere. 

A  natural  framework  for  this  approach  is  the  equilibrium  statistical  me¬ 
chanics,  since  it  deals  with  systems  with  many  degrees  of  freedom  described 
by  Gibbs  probability  distributions.  In  principle,  statistical  mechanics  allows 
us  to  derive  the  mean  statistical  values  of  the  field  (mean  field  equations)  as 
explicit  functions  of  the  data  and  the  parameters  of  the  model:  an  algorithm 
implementing  this  would  consist  just  of  one  step.  Since  the  analytical  com¬ 
putations  required  to  obtain  these  explicit,  expressions  are  usually  loo  hard, 
some  approximations  have  to  be  made.  We  use  the  mean  field  approximation, 
a  well  known  statistical  mechanics  tool,  to  obtain  an  approximated  solution, 
that  is  given  in  implicit  form  by  a  set  of  non  linear  equations.  We  call  these 
equation  deterministic  to  underline  the  deterministic  character  of  the  whole 
procedure. 

We  concentrate  eur  attention  on  the  partition  function  Z,  th  t  is  the  sum  of 
the  probability  distribution  over  all  the  possible  field  configurations,  since  it  is 
known  to  contain  all  the  information  about  the  system.  The  idea  underlying 
our  approach  is  first  to  eliminate  the  line  process  degrees  of  freedom  from  Z: 
we  will  show  that  in  so  doing  the  effect  of  their  interaction  with  the  field  can 
be  simulated  by  a  temperature  dependent  “effective  potential”  that  depends 
only  upon  /.  Its  use  is  fundamental  in  the  derivation  of  the  deterministic 
equations,  and  gives  useful  insights  on  the  role  and  the  significance  of  the 
parameters. 

An  advantage  of  such  an  approach  is  that  the  solution  of  the  deterministic 
equations  is  faster  than  the  Monte  Carlo  techniques,  fully  parallelizable  and 
feasible  of  implementation  on  analog  networks.  The  possibility  of  writing  a 
set  of  equations  is  also  useful  for  a  better  understanding  of  the  nature  of  the 
solution  and  of  the  parameters  of  the  model. 

We  discuss  two  different  MRFs  models.  The  energy  function  of  the  first 
model  has  been  already  studied  by  several  authors[2][l l]i  13](!2j.  It  is  inter¬ 
esting  to  notice  that  the  GNO  algorithm,  proposed  by  Blake  and  Zisserman 
[2]  in  ad  hoc  manner,  arises  naturally  in  the  framework  of  statist!-  a]  mechan¬ 
ics.  This  estabilish  a  connection  between  MRFs  and  determinist  ic  algorithm 
already  used  in  vision. 

In  the  second  model  we  define  an  energy  function  with  an  extra  term,  that 
establishes  an  interaction  between  the  line  process  at  neighborhood  sites. 
This  interaction  can  stimulate  the  creation  of  a  line  at  a  particular  site  if 
a  line  at  a  neighborhood  site  has  been  created.  This  term  allows  the  data 

■> 


to  be  smoothed  and.  at  the  same  time,  the  contrast  at  the  discontinuities 
to  be  enhanced.  Typical  edge  detection  features  like  hysteresis,  threshold 
and  suprathreshold,  which  do  not  arise  from  the  previous  energy  function, 
will  arise  naturally  from  the  model.  We  point  out  that  this  model  is  a  gen¬ 
eralization  of  the  previous  one,  that  can  be  recovered  by  simply  setting  an 
appropriate  parameter  to  zero. 

These  two  methods  can  be  applied  to  dense  data  and,  with  small  modi¬ 
fication,  to  sparse  data  as  well.  The  problem  of  surface  reconstruction  (and 
image  restoration)  from  sparse  data  is  addressed  and  an  algorithm  to  perform 
these  tasks  is  obtained  and  implemented.  We  also  outline  an  algorithm  that 
solves  the  problem  of  aligning  the  discontinuities  of  different  visual  models 
with  intensity  edges  via  integration. 

The  paper  is  organized  in  the  following  way:  Chapter  2  presents  an  overview 
of  MRFs  in  vision.  Chapter  3  discusses  the  deterministic  approximation  of 
MRFs  for  the  three  energy  functions  mentioned  above.  Chapter  4  shows 
how  to  estimate  the  parameters  of  the  models.  In  Chapter  5  some  results  are 
described.  Chapter  t  discusses  applications  for  sparse  data  and  integration 
of  visual  modules  with  intensity.  Chapter  7  concludes  the  paper. 


2  MRFs  for  smoothing  fields  and  detecting 
discontinuities 

Here  we  briefly  summarize  how  MRFs  have  been  used  in  computer  vision.  A 
more  extensive  discussion  is  given  in  Geman  and  Geman  [9],  Marroquin[l2], 
Chou[4],  Gamble  and  Poggio[6],  Gamble.  Geiger,  Poggio,  Weinshall  [5]. 
Consider  the  problem  of  approximating  a  surface  given  sparse  and  noisy 
depth  data,  on  a  regular  2D  lattice  of  sites.  We  think  the  surface  as  a  field 
(surface-field)  defined  in  the  regular  lattice,  such  that  the  value  of  this  field 
at  each  site  of  the  lattice  is  given  by  the  surface  height  at  this  site.  The 
Markov  property  asserts  that  the  probability  of  a  certain  value  of  the  field  at 
any  given  site  in  the  lattice  depends  only  upon  neighboring  sites.  According 
to  the  Clifford-Hammersley  theorem,  the  prior  probability  of  a  state  of  the 
field  /  has  the  Gibbs  form: 


P(f)  =  — 

*1 


(2.1) 


where  /  is  the  field,  e.g.  the  surface-field,  Zj  is  the  partition  function,  U(f)  = 
•£'«(/)  is  an  energy  function  that  can  be  computed  as  the  sum  of  local 
contributions  from  each  lattice  site  and  f3  is  a  parameter  that  is  called  the 
inverse  of  the  natural  temperature  of  the  field.  If  a  sparse  observation  g  for 
any  given  surface-field  /  is  given  and  a  model  of  the  noise  is  available  then 
one  knows  the  conditional  probability  P(g\f).  Bayes  theorem  then  allows  to 
write  the  posterior  distribution: 


P(f  \g)  = 


P[(l\f)P{.t )  _  J_  -DEU\a) 

P(g)  Z 


(2.2) 


As  a  simple  example,  when  the  surfaces  (surface-fields)  are  expected  to  be 
smooth  and  the  noise  is  Gaussian,  the  energy  is  given  by 


E(f\g)  =  -  9i?  +  “£(/,  -  f,?\,  (2.3) 

>  JtN, 

where  7 ,•  =  1  or  0  depending  on  whether  data  are  available  or  not  and  ,V,  is 
a  set  of  site*  in  an  arbitrary  neighborhood  of  the  site  i.  The  maximum  of 
the  posterior  distribution  or  other  related  estimates  of  the  “true”  data-field 
value  can  not  be  computed  analytically,  but  sample  distributions  of  the  field 
with  the  probability  distribution  of  (2.2)  can  be  obtained  using  Monte  Carlo 
techniques  such  as  the  Metropolis  algorithm  [14].  These  algorithms  sample 
the  space  of  possible  values  of  the  surface-field  according  to  the  probability 
distribution  P(f\g)- 

One  of  the  main  attractions  of  MRFs  models  is  that  they  can  deal  directly 
with  discontinuities.  Geman  and  Geman  [9]  introduced  the  idea  of  another 
field,  the  line  process,  located  on  the  dual  lattice,  and  representing  explic¬ 
itly  the  presence  or  absence  of  discontinuities  that  break  the  smoothness 
assumption  (2.2).  The  dual  lattice  is  another  lattice  coupled  with  the  data 
field  lattice  such  that  for  each  site  of  the  data  field  lattice  there  are  two  sites 
of  the  dual  lattice,  one  site  corresponding  to  the  vertical  line  and  the  other 
one  to  the  horizontal  line.  The  associated  prior  energy  then  becomes: 

Ei(f,  /)=  £  (/,  -  mi  -  lij)  +  £  Vcihj)  (2.5) 

jt  tv.  c 

where  l,,  is  the  element  of  the  binary  held  /  located  between  site  j.  The  term 
Vc ( /,j ) ,  where  C  is  a  clique  defined  by  the  neighborhood  system  of  the  line 
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process  (binary  field),  reflects  the  fact  that  certain  configurations  of  the  line 
process  are  more  likely  to  occur  than  others.  Depth  discontinuities  are  usually 
continuous  and  non-intersecting,  and  rarely  consist  of  isolated  points.  These 
properties  of  physical  discontinuities  can  be  enforced  locally  by  defining  an 
appropriate  set  of  energy  values  Vc{l,j)  for  different  configurations  of  the 
line  process  ([9],  [13], [6].)  In  our  models  the  cliques  will  be  simplified  to  the 
nearest  neighbors. 

Two  basic  problems  have  arisen  in  using  MRFs  to  solve  vision  problems.  The 
first  one  is  the  amount  of  computer  time  used  in  the  Metropolis  algorithm 
or  in  simulated  annealing  [10].  The  second  problem  is  to  how  estimate  the 
parameters  of  the  energy  function  since  they  have  been  estimated  in  an  ad 
hoc  manner. 

We  propose  to  approximate  the  solution  of  the  problem  formulated  in  the 
MRFs  frame  with  its  “average  solution.”  The  mean-field  theory  (which  is 
explained  in  the  next  chapter)  allows  us  to  find  deterministic  equations  ana¬ 
lytically  for  MRFs  whose  solution  approximates  the  solution  of  the  statistical 
problem. 

3  A  deterministic  approximation  of  MRFs 

3.1  The  Line  Process  for  two  dimensions 

The  line  process  was  described  in  chapter  2  for  the  one  dimensional  case. 
Now  we  generalize  it  to  two  dimensions.  In  this  case  we  define  a  horizontal 
line  process  hXJ  and  a  vertical  line  process  i/tJ.  The  line  process  connects 
the  site  (i,j)  to  the  site  (i,j  —  1),  while  Vi:  connects  the  site  ( i,j )  to  the 
site  (i  —  1  ,j).  This  is  illustrated  in  figure  1.  It  is  important  to  notice  that 
the  horizontal  line  process  is  determined  by  the  gradient  of  the  field  in  the 
vertical  direction,  and  the  vertical  line  process  is  determined  by  the  gradient 
of  the  field  in  the  horizontal  direction.  In  this  way  we  can  reduce  the  two 
dimensional  problem  to  two  one  dimensional  problems,  provided  that  the 
horizontal  and  vertical  line  processes  do  not  interact.  At  the  end  we  can  add 
both  line  processes  to  get  the  edge  map.  This  approach  will  be  exploited  in 
the  next  chapters.  We  point  out  that  the  fields  ftJ ,  hXJ  and  are  defined 
in  the  same  lattice  instead  of  having  ft]  in  one  lattice  and  htJ,  i\j  in  a  dual 
lattice. 


Figure  1:  The  surface  field  f,  the  horizontal  line  process  h  and  the  vertical 
line  process  v  are  n  presented  in  two  sifts,  (i,j)  and  (i+..\  of  the  lattice. 

There  are  three  fields  defined  in  each  silt  of  the  lattici  as  opposed  to  harimj 
one  field  in  a  lattici  and  the  othrr  two  fields  m  a  dual  lattici. 
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3.2  The  Weak  Membrane  energy 

Or  how  to  smooth  the  field  but  not  at  the  discontinuities 
The  Weak  Membrane  energy  is  given  by 

Ei(f,  h,v)  =  Ejg{f)  +  £■//(/,  h,v)  +  Efih.v)  (3.1) 

where 


E/g(f)  =  Z(f,,J-gi,})2  (3.1a) 


Eji(f,h,v )  =  ~  fi,j- i)2(l  -  hi,})  +  {/,,}  -  fi- i,j)2(l  -  (3.16) 

E,(h,v)  =  7  Yi(hi,j  +  w,j)  (3.1c) 

and  a  and  7  are  positive  valued  parameters. 

The  first  term,  as  in  the  previous  case,  enforces  closeness  to  the  data 
and  the  second  one  contains  the  interaction  between  the  field  and  the  line 
processes:  if  the  horizontal  or  vertical  gradient  is  very  high  at  site  (i,  j ) 
the  corresponding  line  process  will  be  very  likely  to  be  active  (6, j  =  1  or 
Vi^  =  1),  to  make  energy  decrease  and  signal  a  discontinuity.  The  third  term 
takes  into  account  the  price  we  pay  each  time  we  create  a  discontinuity  and 
is  necessary  to  prevent  the  creation  of  discontinuities  everywhere. 


3.3  Mean  field  theory  and  Weak  Membrane 

We  assume  that  there  is  uncertainty  in  the  model  and  then  the  statistical 
framework  can  be  used  to  incorporate  t l,e  uncertainty  as  described  by  f'2.1  1 
Using  mean  field  theory  it  can  be  shown  (see  Appendix  A)  that 


T~  _  1  c/F 


(-1.2) 


where  F  —  ~^lnZ  is  the  free  energy.  Z  is  the  partition  function  defined  as 


1 


Z  = 

If} 

where  H{/}  means  the  sum  over  all  the  possible  configurations  {/}  of  the 
system. 

The  computation  of  the  partition  function  is  equivalent,  in  this  case,  to  the 
evaluation  of  a  multi  dimensional  integral  which  can  not  be  explicitly  solved, 
due  to  the  interaction  between  all  the  variables.  Even  if  an  exact  solution 
can  not  be  found,  we  can  st.ill  obtain  a  good  solution  making  use  of  the  so 
called  mean  field  approximation.  The  mean -field  approximation  is  a  general 
tool  used  in  statistical  mechanics  that  consists  in  substituting  the  interaction 
among  the  fields  at  different,  locations  by  the  interaction  of  the  iieid  at  each 
site  with  the  mean  field  value  at  different  locations  (see  appendix  A) 

After  this  step  the  partition  function  factors  into  the  product  of  single-site 
partition  functions,  that  can  usually  be  computed 


Z"'! 


nC' 

•  a\ 


-HE" 


Now  the  partition  function,  and  then  the  free  energy  are  functions  of  the 
mean  values  of  the  fields,  which  are  still  unknown.  These  values,  however, 
are  usually  related  to  the  derivatives  of  the  free  energy  with  respect  to  some 
external  field,  as  in  (3.2).  Substituting  the  free  energy  with  its  mean-field 
approximation,  we  obtain  a  set  of  implicit  and  self-consistent  equations  for 
the  fields,  usually  called  mean-field  equations.  In  this  sense  the  mean-field 
equations  give  a  deterministic:  solution  to  the  MRFs  problem,  since  they 
relate  explicitly  the  mean  values  of  the  field  to  the  data.  They  could  be 
implemented  by  a  deterministic  network. 

In  the  ease'  of  (3.1)  the  partition  function  becomes 


(3.3) 


U) 


it  , 


where  O’,  ,  —  '  —  n 


fV;„,  -  '■  -a;./,  a*,  -  /, 


and  A'  ,  — 


j- i- 
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3.3.1  Averaging  Out  the  Line  Process 

The  contribution  of  the  line  process  to  the  partition  function  can  be  exactly 
computed.  Indeed  the  line  process  term  in  (3.3)  is  the  partition  function 
of  two  spin  systems  (h  and  v)  in  an  external  field  ( Gh  and  Gv)  with  no 
interaction  between  neighboring  sites.  Then  each  spin  contributes  to  the 
partition  function  independently  from  the  others  and  its  contribution  is  (1  + 
e~0G''i)  for  the  horizontal  field  and  a  similar  factor  for  the  vertical  one.  The 
partition  function  can  then  be  rewritten  as 

Z  =  +  e-^)(1  +  e-^)  (3.4)  . 

{/}  a 

The  capability  of  computing  the  line  process  contribution  to  Z  allows  us  to 
obtain  a  solution  for  the  mean  values  of  the  fields  k  and  v  once  a  solution  for 
the  field  is  provided.  As  an  intermediate  step,  after  some  algebra  derived  in 
appendix  B  we  obtain 

=<  i77a>;  >  md  =<  i"+  e«:„  >  (3-5) 

where  <  ...  >  means  the  statistical  mean  value.  In  order  to  get  a  deterministic 
solution  for  the  line  process  we  need  to  do  some  approximations.  We  will 
assume  then  that  we  can  neglect  the  statistical  fluctuations  of  the  field  and 
so  we  replace  the  value  of  the  field  in  (3.5)  with  its  mean  value.  This  is  in 
essence  the  mean-field  approximation. 


1  _|_  gPh-adi.j-f, 


Vij  = 


1  +  e0h-a(L,}-I,,,-xY) 


(3.6). 


In  the  zero  temperature  limit  (ft  — »  oo)  (3.6)  becomes  the  Heaviside  function 
(1  or  0)  and  the  interpretation  is  simple:  when  the  horizontal  or  vertical 
gradient  (f{j  -  or  are  larger  than  a  threshold  (^)  a 

vertical  or  horizontal  discontinuity  is  created,  since  the  price  to  smooth  the 
function  at  that  site  is  too  high.  This  leads  to  a  clear  interpretation  of  the 
parameter  7,  as  it  will  be  discussed  in  section  4.2. 

We  notice  that  (3.6)  can  account  for  diagonal  lines.  We  show  with  the 
example  of  a  rectangle  in  a  square  lattice  (see  figure  2a).  The  threshold 


diagonal  linos 
horizontal  +  vortical  linos 


Figure  2:  a)  At  the  discontinuity,  A^2  =  (B  —  A)2  and  A,”,-2  =  0.  A 
chain  of  horizontal  line  processes  are  created  for  (^  <  (B  —  A)2),  b)  At  the 
discontinuity,  A^2  =  (B  —  A)2  and  AJj2  ~  [B  —  A)2 .  A  chain  of  horizontal 
and  vertical  line  processes  are  created  for  (^  <  (B  —  A)2). 

to  create  a  discontinuity  line  is  invariant  under  rotations  of  the  rectangle  with 
respect  to  the  lattice  (see  figure  2b).  The  diagonal  line  is  perceived  from  a 
chain  of  horizontal  and  vertical  lines  (like  a  stair  case).  The  important  result 
obtained  from  (3.6)  is  that  although  the  total  line  created  in  figure  2b  is  y/2 
times  bigger  than  in  figure  2a  the  threshold  to  create  a  line  has  been  kept 
the  same. 

3.3.2  The  Effective  Potential 

We  discuss  how  the  interaction  of  the  field  f  with  itself  lias  changed  after 
the  line  process  has  been  eliminated  from  the  partition  function.  From  (3.4) 
we  notice  that  the  partition  function  can  be  rewritten  as 

(/} 
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Figure  3:  The  effective  potential  is  shown  as  a  function  of  A^.  (For  illustra¬ 
tion  we  set  A "•  =  0.)  a)  For  0  —  0.002.  b)  Zero  temperature  limit  (/3  — ►  oo). 

where 

=  I >(Af/  +  A»/)  -  I/n[(l  +  e'^Kl  +  e'^)] 
and  Ejg(f)  is  given  by  (3.1a). 

This  is  the  partition  function  of  a  system  composed  of  one  continuous  valued 
field,  whose  energy  is  Efg  +  EeJJ-  We  interpret  this  result  as  the  effect 
of  the  interaction  of  the  line  processes  with  the  field  /.  This  effect  can 
be  simulated  by  modifying  appropriately  the  interaction  of  the  field  with 
itself,  substitut'ng  the  smoothing  term  in  the  energy  function  with  a  new 
temperature  dependent  potential. 

In  figure  3  the  effective  potential  is  depicted  for  different  temperatures.  It 
simulates  the  effect  of  the  line  processes  on  the  field  /.  Notice  that  the  energy 
function  is  still  the  sum  of  local  interactions  between  first  neighbors.  For  the 
zero  temperature  limit  one  can  see  in  figure  3  that  the  smoothing  term  is 
active  only  when  the  gradient  is  smaller  than  a  threshold,  proportional  to 
the  ratio  between  7  and  a.  When  the  temperature  is  different  from  zero  the 
border  (threshold)  of  the  smoothing  region  is  no  longer  well  defined  due  to 
thermal  noise.  It  corresponds  to  an  adaptive  threshold  that  depends  on  the 
temperature. 


3.3.3  A  Deterministic  Solution  for  / 

The  introduction  of  the  effective  potential  allows  us  to  obtain  a  set  of  de¬ 
terministic  equations  for  the  field  /.  The  mean  field  theory,  as  explained 
in  section  3.1,  can  be  applied  to  the  partition  function  of  equation  (3.4)  to 
obtain  a  set  of  mean-field  equations  depending  on  the  temperature.  In  terms 
of  statistical  mechanics  the  data  term  plus  the  effective  potential  represent 
the  free  energy  of  the  system.  The  mean  field  solutions  are  obtained  by  min¬ 
imizing  the  free  energy.  The  set  of  deterministic  equations  can  be  written 
as 

*?-(£/,  +  £.//(/»  =  0 
°Ji,j 

and  after  some  computation 

A j  =  9i,j  -  -  /i.j-i )( 1  -  )  +_«(/. j+ 1_~  Aj )( 1  -  t>o+i ) 

-  A,;)  +  tt(/i+u  -  (3.7) 

where  htiJ  and  are  given  by  formula  3.6. 

Equation  (3.7)  gives  the  field  at  site  ?,j  as  the  sum  of  data  at  the  same 
site,  plus  an  average  of  the  field  at  its  neighbor  sites.  This  average  takes  in 
account  the  difference  between  the  neighbors.  The  larger  is  the  difference, 
the  smaller  is  the  contribution  to  the  average.  This  is  captured  by  the  term 
(1  where  l,j  is  the  line  process.  At  the  zero  temperature  limit  (/3  — >  oo) 

the  line  process  becomes  1  or  0  and  then  only  terms  smaller  than  a  threshold 
must  be  taken  in  account  for  the  average.  This  interpretation  helps  us  in 
understanding  the  role  of  the  a  and  7  parameters,  as  it  will  be  discussed  in 
chapter  4.  Notice  that  the  form  of  (3.7)  is  suitable  for  the  application  of  a 
fast,  parallel  and  iterative  scheme  of  solution. 

3.3.4  The  Effective  Potential  and  the  Graduate  Non  Convexity 
Algorithm 

We  have  to  point  out  ♦hat  this  energy  function  has  been  studied  by  Blake 
and  Zisserman  [2],  in  the  context  of  edge  detection  and  surface  interpolation. 
They  do  not  derive  the  results  from  the  MRFs  formulation  but  they  simply 
minimize  the  Weak  Membrane  energy  function,  f  rom  a  statistical  mechanics 
point  of  view  the  mean-field  solution  does  not  minimize  the  energy  function. 
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but  this  becomes  true  in  the  zero  temperature  limit,  so  their  approach  must 
be  recovered  from  the  MRFs  formulation  in  this  limit.  This  is  indeed  the 
case,  and  it  is  easy  to  show  that  the  effective  potential  becomes  the  Blake 
and  Zisserman  potential  when  /3  goes  to  infinity.  In  order  to  obtain  the  min¬ 
imum  of  the  energy  function  E\  Blake  and  Zisserman  introduce  the  GNC 
(Graduate  Non  Convexity)  algorithm,  which  is  different  from  our  determin¬ 
istic  scheme,  but  can  be  embedded  in  the  MRFs  framework  in  a  natural  way. 
Let  us  review  briefly  the  GNC  algorithm.  The  main  problem  with  the  Weak 
Membrane  Energy  is  that  is  not  a  convex  function  and  a  gradient  descent 
method  can  not  be  applied  to  obtain  the  minimum  because  one  could  be 
trapped  in  a  local  minimum.  In  order  to  solve  this  problem  Blake  and  Zis¬ 
serman  introduce  a  family  of  energy  functions  E^v\  depending  continuously 
on  a  parameter  p ,  pe[0, 1],  such  that  E ^  is  convex,  E ^  =  E1  and  E ^  are 
non  convex  for  pc[0, 1).  Gradient  descent  is  successively  applied  to  the  energy 
function  E^  for  a  prescribed  decreasing  sequence  of  values  of  p  starting  from 
p  =  1,  and  this  procedure  is  proved  to  converge  for  a  class  of  given  data.  The 
construction  of  the  family  of  energy  functions  E^is  ad  hoc  and  uses  piece- 
wise  polynomials.  In  our  framework,  a  family  of  energy  functions  with  such 
properties  is  naturally  given  by  Eej/T '  where  T  is  the  temperature  of  the 
system.  The  GNC  algorithm  can  then  be  interpreted  as  the  tracking  of  the 
minimum  of  the  energy  function  as  the  temperature  is  lowered  to  zero  (like 
a  deterministic  annealing).  In  this  way  the  approach  of  Blake  and  Zisserman 
can  be  viewed  as  a  deterministic  solution  of  the  MRFs  problem,  even  if  it 
does  not  fully  exploit  the  possibility  of  obtaining  deterministic  equations  for 
the  surface  and  discontinuity  fields.  The  results  obtained  by  applying  this 
method  to  edge  detection,  for  example,  are  good,  and  a  pattern  of  meaningful 
discontinuities  can  usually  be  recovered  [2].  However,  sometimes  the  full  set 
of  discontinuities  is  not  obtained:  when  the  gradient  of  the  image  brightness 
is  under  the  threshold,  a  discontinuity  may  not  be  detected,  even  though  it 
would  be  necessary,  for  example,  to  close  a  contour. 

If  interaction  between  lines  (self  interaction  of  the  line-field)  is  introduced  in 
the  energy  function  this  problem  can  be  overcome,  as  we  discuss  next. 

Figure  4  summarizes  the  procedure  used  before  to  obtain  deterministic 
algorithms  from  the  statistical  methods. 
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Sketch  of  the  Method 


More  precisely  , 


Compute  Z  and  Mean  Held  equations 

1.  Sum  P  over  1  3.  Sum  over  f,  but  It  Is  bard 

Obtain  Efj(f)  and  So,  find  f  that  minimize  E^(f) 

obtain  mean  field  equation  for  1  (mean  field  approximation) 

Aditlonal  feature:  Deterministic  annealing  (Increase  b) 


Figure  4: 


i 


Overvie  w  of  the  method 
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3.4  Improving  the  weak  membrane  model 

Or  how  to  smooth  the  field  while  enhancing  the  differences  at  the  discontinu¬ 
ities 

So  far  we  have  not  exploited  an  important  physical  constraint  of  images, 
namely  the  smoothness  of  the  discontinuity  field.  Isolated  discontinuities  are 
very  unlikely  to  occur  and,  on  the  contrary,  the  presence  of  a  discontinuity  at 
a  site  makes  more  likely  the  presence  of  a  discontinuity  at  a  neighboring  site. 
This  smoothness  constraint  on  the  discontinuity  field  can  be  incorporated  in 
the  model  by  simply  adding  a  new  term  to  the  energy  function.  Thus,  the 
total  energy  becomes 

E2  =  Ei  +  En 

where  E\  is  given  by  (3.1)  and  we  define  the  new  term 

Eu  =  —tf 

and  e  is  a  new  parameter  whose  exact  meaning  and  estimation  will  be  ex¬ 
plained  in  the  next  section.  To  make  more  evident  the  meaning  of  the  new 
term  we  notice  that 


El  +  E»  =  El(l +  *11  +  «*(*«  - Ay-l)’  +  Ki  - ».-u)2l  (3.8) 

ij 

where  Ei  is  given  by  (3.1c).  From  (3.8)  it  is  clear  that  e  is  related  to  the  degree 
of  smoothness  of  the  discontinuity  field  and  so  must  be  a  positive  number. 
In  order  to  keep  positive  the  price  we  pay  for  creating  a  discontinuity  and 
to  prevent  the  line  processes  from  being  active  everywhere,  e  has  to  be  less 
than  1. 

We  notice  here  that  this  model  is  a  simplified  version  of  other  possible  po¬ 
tentials.  In  particular,  the  neighbor  size  considered  here  is  at  most  of  two 
pixels  where  Gamble  &  Poggio  (1987)  for  example,  have  discussed  more  so¬ 
phisticated  cliques  composed  of  larger  neighbors. 
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3.4.1  Averaging  Out  the  Line  Process 

Again,  as  in  the  Weak  Membrane  energy  case  we  calculate  the  contribution 
of  the  line  process  in  the  partition  function.  The  partition  function  can  be 
written  as 


2=  £{/}e 


r^v 


(3.9) 


where  the  G s  have  been  defined  in  section  3.3.  We  rewrite  the  line  process 
term  by  substituting  by  A,-,,- (£?£_,•  -  Notice 

that  this  change  is  just  a  rearrangement  of  the  terms  on  the  sum  over  i,j  and 
therefore  does  not  change  the  value  of  Z.  We  also  notice  that  the  line  process 
term  of  (3.9)  is  the  partition  function  of  many  one-dimensional  Ising  spin 
models  in  an  external  field.  For  a  constant  external  field  an  exact  solution 
can  be  obtained  by  the  method  of  the  transfer  matrix  (see  Parisi  [15]).  This 
method  can  be  still  applied  to  our  case  (not  constant  external  field),  giving 
a  compact  expression  for  the  partition  function.  However,  the  expression  we 
obtain  is  highly  non  local  for  the  field  /  and  not  useful  for  fast  computations. 
Because  of  it  we  apply  the  mean  field  approximation  and  replace  x  by 
h,j- i,  h,j+i  by  hi,j+\  ,  by  and  vi+1J  by  v,+J<J  in  (3.9).  Notice 

that  at  this  step  the  mean  field  approximations  refer  to  the  line  process  and 
not  to  the  field  /  (that  is  kept  “frozen”).  The  sum  over  the  line  process 
configurations  (3.9)  can  then  be  computed,  giving 


gmf-hv  _  £  e-052t,i[G>.)-9>,if+a(±*}2+&lJ2)] , 


n«(l  +  +  e-<9(cr,,-^--i,;'r’tLl: 


The  partition  function  is  now  similar  to  the  previous  one  (3.4)  and  the  zero- 
temperature  mean-field  equations  for  the  line  process  field  can  then  be  de¬ 
rived  in  the  same  approximation: 

hi  =  )2  -  1  +  +*■■'+■)  and 

• \.i  =  "iMf..,  -  -  1  +  «■>— -’«.f  (3.10) 


ill 


where  crp(x)  =  7^737  is  the  sigmoid  function.  Notice  that  now  equations 
(3.10)  form  a  set  of  non  linear  equations  and  the  solution  for  the  line  process 
is  not  as  simple  as  before. 

3.4.2  Hysteresis,  Suprathreshold  and  Threshold 

The  solution  obtained  in  (3.10)  for  the  line  process  can  deal  with  the  prob¬ 
lem  of  “streaking”.  Streaking  is  the  breaking  up  of  an  edge  contour  caused 
by  fluctuations  above  and  below  a  threshold  along  a  contour  [3].  This  is  a 
common  problem  with  thresholded  detectors,  in  particular  with  the  scheme 
implemented  by  the  Weak  Membrane  energy.  For  the  third  energy  the  thresh¬ 
olding  is  done  with  hysteresis.  This  is  the  way  that  Canny’s  detector  works 
[3].  The  hysteresis  phenomenon  is  evident  in  (3.10).  With  the  creation  of  a 
horizontal  line  at  a  neighbor  site,  say  {i,j  —  1}  (A,,j_i  =  1),  the  energy  neces¬ 
sary  for  creating  a  line  at  a  site  {«,  j }  decreases  by  r?.  On  the  other  hand,  if 
two  lines  are  created  at  {i,j  —  1}  and  {i,j  +  1}  then  the  energy  decreases  by 
ey.  A  low  threshold  (threshold)  and  a  high  threshold  (supratlueshokl)  arise 
naturally.  The  suprathreshold  for  creating  a  line  is  given  by  ^  <  (A-*7)2,  in 
this  case  a  line  is  always  created.  In  the  same  way  the  lowest  threshold  is 
given  by  ^-(1  —  e)  <  (A-^)2  <  in  this  case  a  horizontal  line  will  be  created 
at  site  {i,j}  if  a  horizontal  line  has  been  created  in  the  sites  {i,j  —  1}  and 
{0  +  !}- 

At  higher  temperatures  the  threshold  and  suprathreshold  are  not  so  well  - 
defined,  and  what  we  have  are  adaptive  thresholds  that  depend  on  the  values 
of  the  field.  These  points  will  become  clearer  when  we  discuss  the  effective 
potential. 

We  point  ou^  that  some  hysteresis  may  occur  on  the  Weak  Membrane  en¬ 
ergy.  However  it  will  happen  just  on  special  cases  and  without  much  of 
control.  A  consistent  framework  to  enhance  and  smooth  images  requires  at 
least  two  thresholds.  The  Weak  Membrane  energy  is  a  one-threshold  model 
and  therefore  can  not  fully  exhibit  hysteresis. 

3.4.3  The  Effective  Potential 

By  analogy  to  section  3.3,  the  effective  potential  is  given  by 
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E,n({)=  £j,){a(A‘/  +  Av/) 

-JM(1  (3.11) 


An  exact  computation  by  means  of  the  transfer  matrix  method  shows  that 
this  potential  is  highly  non  local  and  its  practical  use  is  limited.  What  we 
will  do  next  is  to  use  approximation  techniques  to  obtain  a  potential  that  is 
local  ,  up  to  the  first  neighbors. 

Without  losing  any  generality  we  now  dismiss  the  contribution  of  ihe  vertical 
process  since  it  is  independent  of  the  horizontal  process  contribution.  In 
order  to  obtain  an  expression  for  Ef/f  in  terms  of  the  gradient  field  A^  we 
investigate  (3.10).  Since  the  discontinuity  field  is  smooth  we  first  approximate 
the  term  by  and  (3.10)  at  the  zero  temperature  becomes 


Ki  =  0(Q(/«,J  -  k j-i)2  -  7  +  n hi)  (3.12) 

The  solution  of  (3.12)  depends  on  the  value  of  A^2  =  ( )2  as  follows 

'  0  if  (a£)2  <  J(1  -  e) 

hi=i  0  or  l  if  J(l-e)  <(£?)’<  2 

1  if  (AD!>? 


These  results  are  illustrated  in  figure  5.  We  notice  that  in  the  region  ^(1  — 
e)  <  (A*  )2  <  ^  both  solutions  are  possible  ,  thus  reflecting  a  weakness  in 
the  assumption  we  made.  A  possible  way  to  obtain  a  unique  solution  is  to 
consider  an  average  of  the  two  solutions  such  that  the  higher  the  gradient  of 
/,  the  higher  is  the  value  of  h  (more  likely  for  creating  a  line).  We  interpret 
the  values  of  h  between  1  and  0  as  the  likelihood  of  creating  a  line.  We  ap¬ 
proximate  the  solution  of  the  mean  field  equation  for  h,j  at  any  temperature 
(/?)  by  the  analytical  expression 


k,  =  -  ~ln(  1  +  +  i-f„(  1  +  c-*) 


(3.13) 


IS 


Figure  5:  The  three  figures  illustrate  the  solution  of  (3.12)  according  to  dif¬ 
ferent  values  of  the  gradient  field  f  (A £•)  a)The  line  process  is  on  (hij  =  1 ) 
b)  Two  solution  (three  solution  for  temperature  different  then  zero )  c)  Shows 
the  solution  of  no  line-process  ( hij  =  0) 


Figure  6:  hij  as  a  function  of  the  gradient,  field  (A|y).  Two  solutions  of 
(3.13).  a)  0=1  b)  =  0.002 
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Figure  7:  The  third  energy  potential  for  ji  =  0.002  and  different  values  for  the 
parameter  e.  a)  For  e  —  0.9  b)  For  e  =  0  the  potential  becomes  the  effective 
potential  of  the  Weak  Membrane  energy. 

that  has  the  desired  property  and  should  be  a  reasonable  guess  for  the  de¬ 
pendence  of  hij  on  the  gradient  of  the  field  /. 

Figure  6  shows  two  solutions  of  (3.13)  for  different  temperatures.  Substitut¬ 
ing  the  value  of  given  by  (3.13)  in  (3.11)  we  then  obtain  the  third  energy 
potential,  that  is  plotted  in  Figure  7  for  some  values  of  /?  and  e. 

The  effect  of  the  third  energy  potential  on  the  mean  field  /  can  be  understood 
by  analyzing  the  “force”  (derivative)  of  the  potential  (see  Figure  8).  The 
third  energy  effect  can  be  described  as  follows. 

1.  Smooth  the  field  for  small  values  of  the  gradient  ((A-^  ?  <  (A0)2). 
(Positive  derivative  of  the  potential.) 

2.  Enhance  the  discontinuity  for  (A0)2  <  (A,^)2  <  in  other  words, 
forcing  the  value  of  the  field  at  a  site  to  be  pulled  away  from  the  value 
of  the  field  at  the  considered  neighbors.  (Negative  derivative  of  the 
potential.) 

3.  Neither  smoothing  nor  enhancing  the  field  for  ( A^ ) 2  >  -T  (The 
derivative  of  the  potential  is  zero.) 

We  point  out  that  t  he  third  energy  potential  derived  above  is  a  local  approx¬ 
imation  of  the  mean  field  approximation  over  the  discontinuity  field. 
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Figure  8:  The  third  energy  force  for  3=  0.002  and  different  values  of  the 
parameter  t  a) For  e=0.9  .  b)  For  t  —  0  the  force  is  the  same  as  for  the  Weak 
Membrane  energy. 

3.4.4  A  deterministic  solution  for  / 

Here  again  the  same  arguments  of  section  3.3.3  apply  to  obtain  a  determin¬ 
istic  solution  for  the  field  /.  After  some  computations  we  then  obtain: 

fij  =  9i,i  -  *>. j )  +  aAj.j+A 1  ~  l ) 

-otA^l  -  hjj)  +  a&hijil  ~  h+ij) 

+  aeA^hijerpii  -  a(A^Q _ 

-  -Q(^t-n)2) 

+  a(Ai,vm°0{'l  ~  af(AJj-Q _ 

-  a*Avj+1Vij+1<T0(7  -  q(AJ;+1)2)  (3.14) 

4  Parameters 

The  parameters  a,  7,  e  and  ft  must  be  estimated  in  order  to  develop  an 
algorithm  that  smoothes,  enhances  and  finds  the  discontinuities  of  the  given 
data-field. 


4.1  The  parameter  a 

The  parameter  a  controls  the  balance  between  the  “trust’'  in  the  data  and 
the  smoothing  term.  The  noisier  are  the  data  the  less  you  want  to  “trust"  it 
so  a  is  larger,  the  less  noisy  are  the  data  the  more  you  “trust”  it  so  a  should 
be  smaller.  To  estimate  a  various  mathematical  methods  are  available.  The 
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generalized  cross  validation  method  introduced  by  Wahba  [17]  and  the  stan¬ 
dard  regularization  method  described  by  Tikhonov  [1G,  1]  were  studied  by 
Geiger  and  Poggio  [7]  to  estimate  a.  The  standard  regularization  method 
gives  [7]  an  equation  relating  a  with  the  signal-to-noise  ratio  that  for  the 
continuous  case  is 


rr 

Jo  Jo 


SU\V)( u.-?  +  »z‘)2  p  p  A 

m  J0  p 


A ;(^n)(vC‘  -t  i'2 ) 


.2  ,,2  \  13 


( 1  +  rrfu.-2  4-  r2)) 


-rL :  du  (4.1) 


where  is  the  spectral  density  (power  spectrum)  of  the  signal  f(x.y) 

and  N(u>,  v)  is  the  spectral  density  of  the  noise  v(r,y),  assuming  that  v(x,y) 
is  stationary. 

When  an  estimation  of  the  noise  is  not  available,  it  is  necessary  to  use 
the  generalized  cross  validation  method  (GCY).  GCV  states  that  the  opti¬ 
mal  value  of  q  can  be  obtained  by  minimizing  the  functional  (here  in  one 
dimension) 

,•/  1  p  [./'V  :>(h)  _  </.]*  ,  /  .  .-o 

l  (c»)  =  -2^77 - 1 — k[C  ^  l4-2) 

n  7TTi  (1  -  akk{oc)Y 

where  fn_a(tt)  is  the  smoothed  solution  given  by  (3.1). 


do ))/( 1 - 5Z  <1jj^ 


(O;  U-  "u 


and  nkk( o)  =  )(h ! 

For  the  purpose  of  smoothing  images  both  methods  give  about  the  same 
values  of  a  [7].  Better  results  should  be  obtained  if  the  estimation  of  o  is 
done  locally  (for  a  small  neighborhood). 

4.2  The  parameter  p 

From  (3.6)  one  can  see  that  ,p  is  tin*  threshold  for  creating  a  line  in  the 

Weak  Membrane  energy,  bet’s  cal!  £  the  value  .  From  1 1n*  expression  of 
the  effective  potential  we  notice  that  i!  the  gradient  A f:,  is  above  £  there  is 
no  smoothing  and  i!  the  gradient  is  Ix  low  £  then  smoothing  is  applied.  The 
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parameter  £  defines  the  resolution  of  the  system  and  we  explain  this  with 
two  examples. 

In  the  first  example  suppose  we  are  working  with  the  stereo  module,  so 
the  data  field  is  a  depth-field.  In  this  case  £  is  the  threshold  for  the  changes 
in  depth  to  be  called  a  depth  discontinuity.  This  value  is  determined  ac¬ 
cording  to  the  resolution  of  the  stereo  system  available  and/or  to  the  desired 
resolution  that  one  is  interested. 

In  the  second  example  the  field  is  the  intensity  data  and  the  parameter  f 
represents  the  threshold  for  detecting  edges.  This  value  is  somehow  arbitrary, 
and  probably  context  dependent.  A  variation  of  £  =  16  units  to  £  =  32 
units  for  an  8-bit  array  image  is  likely  to  be  in  agreement  with  the  human 
threshold.  The  exact  value  of  £  depends  on  the  attention  of  the  observer 
and/or  the  sensitivity  of  the  system. 

For  the  second  energy  the  absolute  threshold  to  create  an  edge  is  also 
given  by  £.  However  when  there  is  support  from  neighbor  edges  this  threshold 
can  be  lowered  to  £  x  (w( T  —  e)  (see  (3.10)).  This  suggests  setting  7  so  to 
guarantee  that  the  highest  noise  gap  is  smoothed  (we  are  assuming  that  noise 
gaps  are  isolated  features).  In  this  case  the  e  parameter  will  be  set  to  assure 
that  at  the  edges  the  images  will  not  be  smoothed  but  perhaps  enhanced.  In 
this  case  the  value  of  £  =  30  may  be  desired  for  an  8-bit  array. 

4.3  The  parameter  € 

The  parameter  e  allows  the  energy  to  be  more  general  by  controlling  the 
amount  of  propagation  of  the  line.  So,  once  a  line  is  created,  the  price  to  pay 
for  creating  another  line  next  to  it  wiii  !,e  lowered  by  the  amount  of  7 e.  In 
other  words,  from  (3.1 1)  one  can  see  that  the  difference  in  the  energy  for  when 
a  line  has  been  created  and  when  no  line  has  been  created  at  pixel  (?  —  l,  j). 
is  given  by  This  is  what  characterizes  the  threshold  and  suprathreshold 
or  the  hysteresis  phenomena  [3].  The  threshold  is  given  by  ar>d  the 

suprathreshold  by  y The  parameter  f  varies  from  0  to  1.0  and  when  is  zero 
reduces  the  third  energy  to  the  Weak  Membrane  energy.  When  e  =  1.0  lines 
are  created  everywhere,  since  once  a  line  is  created  there  is  no  cost  in  creating 
another  one  and  then  it  propagates  indefinitely  (to  the  image  limits).  How 
much  does  one  want  to  propagate  a  line?  How  much  should  the  difference 
between  the  threshold  and  the  suprathreshold  be?  Which  exact  value  of  e 
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should  be  chosen?  According  to  the  discussion  above  about  7  we  conclude 
that  e  should  be  chosen  (assuming  we  already  have  a  and  7)  to  guarantee 

that  below  the  desired  edge  threshold.  Again  it  depends  on  what 

one  wants  to  do  with  the  data.  For  detecting  edges  of  objects  in  a  scene  one 
wants  to  have  high  suprathreshold  and  threshold  (usually  obje  *  boundaries 
exhibit  high  gradients)  and  e  large  (bigger  than  0.5)  so  that  all  the  object 
boundaries  are  detected,  included  the  exceptional  boundary  pixels  with  a 
smaller  gradient. 

4.4  The  parameter  3 

The  parameter  0  controls  the  uncertainty  of  the  model.  The  smaller  is  0 
the  more  inaccurate  is  the  model.  This  suggests  that  for  solving  the  mean 
field  equations  a  rough  solution  can  be  obtained  for  a  small  value  of  0  (high 
uncertainty)  and  thereafter  we  can  increase  0  (small  uncertainty)  to  obtain 
more  accurate  solutions.  This  can  be  called  deterministic  annealing. 

We  also  notice  that  0  multiply  the  first  term  on  (3.1)  to  give 

which  suggests  that  the  inverse  of  0  gives  the  t  tandard  deviation  of  the  noise. 
We  keep  in  mind  that  0  also  multiplies  the  other  parameters  of  the  model 
and  therefore  by  estimating  0  b\  the  amount  of  noise  the  estimation  of  the 
other  parameters  have  to  be  scaled  by  this  factor. 

5  Results 

5.1  Implementation 

For  the  implementation  the  zero  temperature  limit  equations  have  provided 
results  as  good  as  the  deterministic  annealing  with  a  faster  computational 
time.  We  do  not  have  proof  of  convergence  but  only  suggestive  experimental 
results. 

In  order  to  find  the  mean  field  solution  for  the  third  energy,  we  solved  (3.14) 
together  with  (3.10)  in  a  coupled  and  iterative  way.  Move  precisely: 

For  simplicity,  let's  write  (3.14)  as 


fi,j  —  P'ifjt.ji  ft,j  i  '  fi.j—  1  •  fi+l .  M  fi,j  +  \  i  hi,j »  V i,j —  1  ) 

and  (3.10)  as 


H (/>,}  •  ft,}- 1 1  h‘,J- 1 t  +  1  )  l  ‘,J  V(fi,j  >  /i-l  »  vi-i  ,j  1  vi+1  ,j  ) 

The  algorithm  solves  these  equations  iteratively  by  updating  (for  the  parallel 
case)  in  the  following  way 


n;'  -  -  da,, . ft,  («.d 


and 


vr  =  t,Knt)  ttttd  tft'  = 


where  a;  controls  the  rate  of  change  and  the  index  n  indicates  the  step  of  the 
iterative  process. 

For  a  serial  implementation  we  first,  update  the  even  sites  (like  the  white 
squares  in  a  chess  board)  and  then  the  odd  sites  (like  the  black  squares). 
Typically  the  algorithm  has  converged  in  10  iterations  which  takes  about  1 
minute  for  images  of  64  x  64  pixels  on  a  Symbolics  3600. 


5.2  Synthetic  Images 

We  first  tested  the  algorithm  on  a  synthetic  step  edge  image  with  noise  added. 
This  is  a  good  test-image  since  locally  many  edges  on  real  images  are  of  this 
type. 

The  step  edge  image  is  an  8-bit  array  of  64  x  64  pixels  with  a  step  intensity 
of  140  units  (see  Figure  9a).  White  noise,  with  standard  deviation  30,  is 
added  to  the  step  edge  (see  Figure  9b).  We  apply  the  algorithm  described  by 
(5.1)  to  reconstruct  the  original  image.  First  we  use  t  =  0.  in  order  to  have 
the  Weak  Membrane  energy.  We  set  o  =  4  (typical  smoothing  parameter) 
and  7  =  24000  (so  that  we  do  not  smooth  the  edge).  The  result  after  20 
iterations  is  shown  in  Figure  9c.  We  notice  that  the  step  edge  starts  to  be 
smoothed  before  all  the  noise  has  been  smoothed.  In  Figure  9d  we  set  e  to  be 
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0.9  and  7  =  45000,  and  the  result  is  more  striking.  Because  we  set  the  ratio 
~  to  be  higher  than  that,  all  the  noise  has  been  smoothed  away.  However,  at 
the  edge  there  is  support  from  the  neighbor  edges  to  lower  the  threshold  to  a 
factor  (1  —  e)  =  0.1.  Therefore  the  noise  is  suppressed  and  the  edges  are  not 
smoothed  but  on  the  contrary  are  enhanced.  We  conclude  that  for  a  noisy 
step  edge,  the  third  energy  is  able  of  retrieving  and  enhancing  the  edge  (see 
figure  9  Id)  while  the  Weak  Membrane  energy  retrieves  the  edge  with  worse 
performance  (see  figure  9c). 

5.3  Real  Images 

When  we  apply  the  second  energy  algorithm  to  a  real  still  life  image  the 
result  is  an  enhancement  of  specular  edges,  shadow  edges  and  some  other 
contours  while  smoothing  out  the  noise  (see  Figure  10).  This  result  is  ob¬ 
tained  consistentently  in  all  the  images  we  used  and  also  with  other  fields 
besides  intensity,  e.g.  ‘‘color”  images.  In  our  case  “Color”  images  are  8-bit 
images  representing  the  ratio  between  the  red  array  and  the  sum  of  red  and 
green  arrays. 

6  Surface  reconstruction  and  alignment  via 
integration 

6.1  Sparse  Data  and  Surface  Reconstruction 

Until  now  we  dealt  with  dense  data  defined  on  allthe  sites  of  the  lattice. 
This  is  not  always  the  case,  especially  if  the  field  is  the  output  of  a  stereo  or 
motion  module.  In  that  case  data  are  given  on  a  subset  of  the  lattice  sites 
and  the  data-field  interaction  term  (Ejg)  in  the  energy  function  has  to  be 
modified  in  the  following  way 

EJa  —  TM.,  ~  9>,j)  ^  rjf,,  ~~  9*.})  ~u ../• 

'•j  «.j 

Here  7,^  is  a  flag  that  is  one  if  a  datum  is  given  at  site  (i,j)  and  0 
otherwise.  This  slight  modification  lias  no  effect,  on  the  theoretical  results, 
and  only  some  changes  take  place  m  the  deterministic  equations  for  the 
surface  field:  the  term  enforcing  closeness  to  data  disappears  when  a  datum 
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Figure  9:  a)  Step  edge  with  IfO  grey  value  units  far  the  step,  lb)  White 
noise  with  standard  deviation  30  grey  value  units  has  been  added,  c)  The 
noisy  image  afire  20  iterations  for  a  =  i,  7  =  2-! 000.  <  —  0.  d)  The  noisy 
images  after  20  /'<  ration  for  a  =  4.  7  =  -15000,  <  --  U.!1 


Figure  10:  a)  Still  lifi  image  Id S  x  IdS  pix'is.  hi  Tht  image  smoothed  with 
e  —  0.5,  7  =  1400  and  o  —  4  for  9  iterations  c)  The  line  process  field  (without 
thinning). 

is  not  given  at  a  bite.  We  rewrite  here  the  deterministic  solution  for  the  field 
/  in  the  case  of  the  Weak  Membrane  hnergy,  since  the  third  energy  equations 
are  modified  in  the  same  way.  In  the  case  of  sparse  data  (3.7)  becomes  then 

/ij7»j  =  -  <>(/.,,  -  /...,  i  K  1  -  +  c\(f,j  4,  --  f.j)(l  - 

"<v(/t.;  -  /--I  j  H  1  -  I ,)  )  1  oi  /,+  !,;  —  ./,,;)(!  —  h,  +  }  t) 

To  apply  this  algorithm  one  needs  first  to  fill  in  the  data.  We  chose  to  fill 
in  the  data  by  averaging  next  neighbors  and  applying  the  algorithm  at  the 
same  time.  So  at  each  step  of  the  algorithm  each  lattice  site  is  visited:  if 
there  is  no  field  value  the  average  neighbor  value  is  taken;  otherwise  we  apply 
the  above  algorithm.  We  notice  that  no  action  is  taken  il  at  a  particular  site 
there  is  no  field  value  and  no  neighbor  held  value. 

From  one  face  image  we  produced  sparse  data  by  randomly  suppressing  70 
%  of  the  data  (see  Figure  II).  Wo  then  applied  the  third  energy  algorithm 
to  sparse  data.  The  parameters  were  kept  the  same  as  the  other  real  image. 
The  reconstruction  from  sparse  data  can  he  applied  to  depth  data  in  which 
case  it  is  usually  called  surface  reconstruct  ion.  We  used  a  stereo  algorit  hm 
based  on  zero  <  r<  >ssmgs  to  obtain  depth  data  from  the  image  shown  below. 
Then  we  applied  the  algorithm  to  reconst  i  net  the  depth  surface.  The  param¬ 
eters  were  different  ,  err. 'ding  to  t  he  .a  bi'i  I,'  for  1  ptli  discontinuity. 
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Figure  11:  a)  A  face  image  of  8-bits  and  128  X  128  pixels,  b)  Randomly 
chosen  70  %  of  the  original  image.  For  display  the  other  30%  are  filled  with 
white  dots,  c)  The  algorithm  described  above  is  applied  to  smooth  and  fill  in 
at  the  same  time  with  e  =  0.9,  7  =  1400  and  a  —  4  for  70  iterations. 
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Figure  12:  a)  The  left  image  (8-bits  and  128  x  128  pixels),  b)  The  right 
image  of  the  same  scene  c)  Sparse  depth  data  obtained  bg  the  stereo  algorithm 
based  on  zero  crossing.  Intensity  represents  depth  calm  s,  d)  The  algorithm 
is  applied  with  e  =  0.0,  -7  =  400  and  a  =  4  to  reconstruct  the  surface,  e) 
Canny  edges  for  figure  a),  f)  The  alignment  algorithm  is  applied  with  the 
same  parameters  as  d)  and  6  =  084.  Depth  discontinuities  aligned  with  the 
intensity  edges  are  obtained. 


6.2  Alignment  of  visual  modules  with  intensity  edges 

The  integration  of  different  visual  modules  to  improve  the  detection  of  the 
discontinuities  can  also  be  addressed  in  this  scheme.  As  suggested  by  Gamble 

6  Poggio  [6]  ,  we  can  add  the  term  6 (v{j  +  /i,j)(l  —  e,j)  to  the  Weak  Membrane 
Energy  or  the  second  energy.  Here  is  an  external  field,  for  example  the 
edge  map  that  is  coupled  with  the  stereo  field.  For  implementation  purposes 
the  only  consequence  of  adding  this  term  is  the  change  of  the  global  parameter 

7  into  the  local  parameter  7 =  7  —  6(1  —  etJ ) .  R.  Thau  implemented 
this  scheme  in  the  Connection  Machine  using  Canny’s  intensity  edges  for 
e,j.  We  first  took  a  pair  of  images  (Figure  12)  and  used  a  stereo  algorithm 
to  find  depth  at  the  zero  crossing  positions.  We  then  used  the  algorithm 
to  reconstruct  depth  everywhere  and  detect  depth  discontinuities.  As  we 
expected  depth  discontinuities  are  aligned  with  the  intensity  edges. 


7  Conclusion 

We  have  used  statistical  mechanics  tools  to  derive  deterministic  approxi¬ 
mations  of  Markov  random  fields  models.  In  particular  we  have  studied  an 
energy  model  that  is  suitable  for  image  reconstruction  or  any  field  reconstruc¬ 
tion.  The  model  has  been  developed  to  include  the  following  characteristics: 

•  the  surface  field  is  smoothed  when  its  gradient  is  not  too  high, 

•  contrast  will  be  enhanced  where  a  discontinuity  occurs  (if  it  is  not  too 
large  already), 

•  the  discontinuity  field  is  likely  to  be  smooth  (isolated  discontinuities 
are  inhibited), 

•  hysteresis  and  adaptive  multiple  threshold  arise  naturally  from  the 
model. 

•  three  parameters  are  needed  to  specify  the  model, 

•  when  one  of  the  parameters  (e)  is  set  to  zero  the  weak  membrane  energy 
is  recovered. 

•  An  understanding  of  the  role  of  the  parameters  is  possible, 
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•  The  model  can  deal  with  sparse  data  and  alignment  of  the  discontinu¬ 
ities  of  different  modules  with  the  intensity  edges. 

We  have  shown  that  the  deterministic  algorithm  of  GNC  can  be  regarded 
as  an  approximation  of  the  gradient  descent  method  with  a  deterministic 
annealing  schedule  to  solve  the  mean  field  equations.  This  suggests  a  unified 
framework  to  connect  different  methods  used  on  image  segmentation,  restora¬ 
tion  and  surface  reconstruction.  We  will  show  in  another  paper[8]  that  several 
deterministic  algorithms  for  image  segmentation  and  reconstruction  are  ap¬ 
proximations  of  two  methods  to  solve  the  mean  field  equations:  the  gradient 
descent  method  discussed  in  this  paper  and  the  parameter  space  method 
discussed  in  [8].  We  derived  a  deterministic  solution  for  the  mean  values  of 
the  surface  and  discontinuity  fields,  consisting  of  a  system  of  coupled  non¬ 
linear  equations.  An  algorithm  has  been  implemented  to  obtain  a  solution 
for  this  system:  it  is  fully  parallelizable,  iterative  and  recursive,  allowing  effi¬ 
cient  computation.  It  would  be  interesting  to  analyze  other  approximations 
or  extensions  of  this  model.  For  example,  a  model  that  includes  interac¬ 
tion  between  the  horizontal  and  vertical  line  processes  could  be  developed  to 
inhibit  self- intersections  of  the  discontinuity  field.  A  term  like  h,j(  1  —  u,j) 
may  be  sufficient.  An  analysis  of  convergence  of  the  algorithm  would  also  be 
important. 
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Appendix  A 

Given  an  energy  function  of  the  form 


eu)  =  £«/,  -  +  IE  m,  />)]> 

i  jcN, 

where  Ni  is  the  set  of  neighbor  sites  of  i,  V(fi,fj)  is  an  interaction  potential 
and  gi  is  an  external  field  (data),  the  partition  function  is  given  by 


Z  =  '£e-0EM 

f 

where  the  E/  means  a  sum  over  all  the  possible  configurations  of  the  system, 
that  is  a  multidimensional  integral  over  all  the  variables  /,,  in  this  case.  The 
mean  field  equation  for  /*,  is  given  by 


(Al) 


From  this  definition,  and  the  definition  of  the  partition  function  Z  ,  the 
following  equality  is  derived: 


1  dZ 
Z  dgk 


=  -2  (3(fk-gk) 


(A2) 


Defining  as  usual  the  free  energy  F  as  F  =  —jjlnZ  we  can  now  obtain,  from 
(A2) 


,  1  OF 

fi~9k  2  dgt 


(/13) 


In  the  case  of  energy  E  the  mean-field  approximation  consists  in  replacing 
the  fields  at  the  neighbor  sites  with  their  statistical  mean  value  Therefore 
the  mean-field  approximation  Em}  to  the  energy  function  becomes 


Em’u)  =  £<(/.  -  s.)2  +  ie  mJ,)} 

t  jtN, 

This  is  a  good  approximation  of  the  original  energy  when  the  fluctuations  of 
the  field  /  are  small,  which  is  the  case  in  our  experiments. 
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Appendix  B 

In  this  appendix  we  will  derive  an  expression  for  the  mean  value  of  the  h 
and  v  fields  of  a  system  described  by  the  Weak  Membrane  Energy 


Noticing  that 


Ei(f,h,v)  =  +  Eji(f,h,v)  +  E,(h,v). 


ExU,h,v)  =  £,(/)  +  +  «wGw) 


where  the  energy  function  Ex  has  been  defined  in  section  1,  the  partition 
function  can  be  written  as 

Z  =  J2e~pElU)  e~0^lhi')G^+Vi-}G^] . 

{/}  {*.»} 

Let  us  consider  just  the  field  h.  By  definition  the  mean  value  hij  is 

hi  =  i  E  Z  (Bl) 

^  {/}  {h.v} 

Then  we  can  rewrite  (Bl)  in  the  following  way: 


hi'i  ~  o  7  1 


-OBiU). 


Zk,{f) 


where 


ZM)  -  £ 

Zhv  represents  the  contribution  of  the  line  process  to  the  partition  function 
of  the  whole  system  for  each  fixed  configuration  of  the  field  /  and  can  be 
seen  as  the  partition  function  of  the  line  process  system  when  the  field  /  is 
kept  “frozen”.  We  notice  that  the  line  process  system  is  simply  a  classical, 
non-interactive  spin  system  in  an  external  field.  Since  there  is  not  interaction 
between  the  spin  variables,  Zhv(f )  can  be  easily  computed,  giving 

ZU/)  =  n<l+<~0O^(l+e-BO'<). 


We  can  now  define  the  free  energy  of  this  spin  system  by  the  usual  relation 


ZUf)  = 

and  substituting  in  (B2)  we  obtain  the  following  equation 

hi  =  ^Z‘-elEM>+r'-M>ig§rFi,Af)-  (03) 

Before  proceeding  further  we  now  notice  that,  having  defined  Fi„(/),  the 
partition  function  of  the  whole  system  can  now  be  rewritten  as 


Z  _  £e-0[E i(/)+Ffc.(/)l. 

{/> 

This  can  be  regarded  as  the  partition  function  of  a  system,  composed  just  by 
the  field  /,  whose  energy  function  is  Ex{f)  -f  Fhv(f)-  Now,  indicating  with 
<  ...  >  the  statistical  mean  value  over  all  the  possible  configurations  of  this 
system,  (B3)  can  be  rewritten  now  in  the  following  way: 


and  then 


an,(/) 

soi 


> 


hj  =< 


1  +  e*%-> 


{BA) 


A  similar  expression  holds  for  Vij.  The  meaning  of  these  equations  is  the 
following:  if  the  field  /  is  kept  “frozen"  the  mean  value  of  the  horizontal 
line  process  is  simply  dgo^K  but  to  obtain  the  true  mean  value  htJ  we  have 

to  average  ^^4^  over  all  the  possible  configurations  of  the  field  /.  whose 

C.  1,3 

probability  distribution  is  given  by 


J£,(/)+Fhv.(/)] 

PU)  =  — j — 

where  the  term  Fhv{f)  is  due  to  the  interaction  of  the  line  process  with  the 
field  f . 
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